!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Alessio Maggioni                   <a.maggioni87@gmail.com>

!>\brief
!!
!! Turbulent compressible channel flow.
!!
!! \n
!!
!! This setting is suitable for the simulation of a turbulent,
!! three-dimensional channel. The viscosity depends on the temperature
!! as
!! \f{displaymath}{
!!  \nu = \frac{\mu_0}{\rho}\left(\frac{T}{T_0}\right)^\alpha
!! \f}
!! for some constants \f$\mu_0\f$, \f$T_0\f$ and \f$\alpha\f$. We then
!! set \f$\nu_T = \frac{\nu}{Pr}\f$.
!!
!! \section dimset Dimensional Setting
!!
!! The test-case setup is specified in terms of some dimensionless
!! numbers. To this end, let us choose some reference values
!! \f{displaymath}{
!!   \rho_r, \quad L_r, \quad V_r, \quad T_r
!! \f}
!! for density, linear dimension, velocity and temperature,
!! respectively. These quantities now define the new <em>units</em>,
!! indicated here as \f$*_u\f$, of the (usual) fundamental dimensional
!! quantities as
!! <center><table>
!!  <tr> <th>variable</th> <th>unit</th> </tr>
!!  <tr>
!!    <td>length</td> <td>\f$L_u = L_r\f$</td> 
!!  </tr> <tr>
!!    <td>time</td> <td>\f$\displaystyle t_u = \frac{L_r}{V_r}\f$</td> 
!!  </tr> <tr>
!!    <td>mass</td> <td>\f$M_u = L_r^3\rho_r\f$</td> 
!!  </tr> <tr>
!!    <td>temperature</td> <td>\f$T_u = T_r\f$</td> 
!!  </tr>
!! </table></center>
!! We then have the equivalences \f$1\,m = \frac{1}{L_r} L_u\f$ for
!! lengths, \f$1\,s = \frac{V_r}{L_r} t_u\f$ for times, \f$1\,kg =
!! \frac{1}{L_r^3 \rho_r} M_u\f$ for masses and \f$1\,K =
!! \frac{1}{T_r} T_u\f$ for temperatures, where the <em>
!! numerical</em> values \f$\rho_r\f$, \f$L_r\f$, \f$V_r\f$, \f$T_r\f$
!! are given in international units. All the dimensional quantities
!! and constants in the model will be expressed with respect to the
!! \f$*_u\f$ units. For instance, given the value
!! \f$\mu_0=1.716\cdot10^{-5}\,kg/ms\f$, the value used for all the
!! computations will be
!! \f{displaymath}{
!!  \tilde{\mu}_0 = \frac{1.716\cdot10^{-5}}{L_rV_r\rho_r} \,
!!  \frac{M_u}{L_ut_u}.
!! \f}
!! Also, denoting by \f$\tilde{U}\f$ the momentum computed by the
!! code, the corresponding value in the international system is
!! \f{displaymath}{
!!  U = \rho_r V_r \tilde{U}\,\frac{kg}{m^2s}.
!! \f}
!! \note A more general approach would also introduce a different
!! number of fundamental (or primitive) variables; this however would
!! introduce excessive complications.
!! 
!! Now that we have chosen our units (at least formally, since the
!! precise values of \f$\rho_r\f$, \f$L_r\f$, \f$V_r\f$ and \f$T_r\f$
!! have not been given yet), we have to prescribe the data of the test
!! case. We proceed as follows:
!! <ol>
!!  <li> the data of the test case are prescribed so that the <em>flow
!!  regime</em>, as identified by a collection of dimensionless
!!  numbers, is of the desired type;
!!  <li> the reference values \f$\rho_r\f$, \f$L_r\f$, \f$V_r\f$ and
!!  \f$T_r\f$ are chosen in order to ensure that the numerical values
!!  of the computed solution are \f$\rho,L,V,T = \mathcal{O}(1)\f$.
!! </ol>
!!
!! Let us thus introduce the following dimensionless numbers
!! \f{displaymath}{
!!  Re = \frac{\rho V L}{\mu}, \quad
!!  Ma = \frac{V}{\sqrt{\gamma RT}}, \quad
!!  Pr = \frac{c_p\mu}{k},
!! \f}
!! where \f$k=c_p\rho\nu_T\f$ is the thermal conductivity.
!! By virtue of our previous requirement n&deg; 2, we obtain (we omit
!! the tildes and assume that all the quantities are expressed with
!! respect to the \f$*_u\f$ units)
!! \f{displaymath}{
!!  Re = \frac{1}{\mu}, \quad
!!  Ma = \frac{1}{\sqrt{\gamma R}}, \quad
!!  Pr = \frac{c_p\mu}{k}.
!! \f}
!! Taking now \f$T_r=T_0\f$, so that \f$\mu\approx\mu_0\f$ (this is
!! not really restrictive, since the ratio \f$T_r/T_0\f$ could be
!! absorbed in the constant \f$\mu_0\f$ anyway), we can determine
!! some problem coefficients as
!! \f{displaymath}{
!!  \mu_0 = \frac{1}{Re}, \quad
!!  R = \frac{1}{\gamma\, Ma^2}, \quad
!!  \nu_{T_0} = \frac{1}{Re\,Pr}.
!! \f}
!! This, together with the ideal gas assumption \f$c_p=c_v+R\f$,
!! completes the specification of the problem, since then we also have
!! \f{displaymath}{
!!  c_p = \frac{1}{(\gamma-1)Ma^2}, \quad
!!  c_v = \frac{1}{\gamma(\gamma-1)Ma^2}, \quad
!!  \kappa = \frac{\gamma-1}{\gamma}.
!! \f}
!! Notice that \f$\gamma\f$ must be specified; for instance, assuming
!! a biatomic gas, we have \f$\gamma=7/5\f$.
!!
!! \note To run the numerical experiments, we do not need to specify
!! any precise value for \f$\rho_r\f$, \f$L_r\f$, \f$V_r\f$ and
!! \f$T_r\f$. This happens because our results are indeed
!! representative of a whole family of similar solutions. Once a value
!! is specified for these variables, all other problem quantities are
!! determined, which corresponds to selecting one specific solution
!! out of a similarity collection.
!!
!! \section random_init Random Initial Condition
!!
!! The initial condition is prescribed as a logarithmic profile with a
!! "random" perturbation. However, we want a perturbation which allows
!! reproducibility of the numerical results, both for different runs
!! and for different grid partitioning. For this reason, we avoid
!! using a random number generator, and instead use a chaotic process
!! to produce a spatially uncorrelated function of the physical
!! coordinates. More precisely, for a given \f$\xi\in[0\,,\,0.5]\f$,
!! we compute
!! \f{displaymath}{
!!   \eta_0 = \xi, \qquad {\rm
!!   for\,\,i=1,20}:\,\,\,\eta_i=r\,\eta_{i-1}(1-\eta_{i-1})
!! \f}
!! with \f$r=3.999\f$.
!!
!! \section forcing_term Forcing Term
!!
!! The forcing term is defined in order to attain the prescribed mass
!! flow \c bar_mass_flow. Since the correct value can not be
!! determined a priori, we use
!! \f{displaymath}{
!!   \underline{f} = \left( \frac{Re_{\tau}^2}{Re^2} - \alpha_P( M_1 -
!!   \bar{M}_1 ) - \alpha_I
!!   \mathcal{M}_{\mathcal{I}_1}\right)\underline{e}_1
!! \f}
!! where the first constant term provides a reasonable guess for the
!! required forcing and the two terms are the proportional and
!! integral corrections, where \f$\underline{M}\f$ is the mass flow as
!! defined in \c mod_dgcomp_flowstate and
!! \f$\underline{\mathcal{M}}_{\mathcal{I}}\f$ is the integrated mass
!! flow as described in the documentation of \c t_dgcomp_stv. The
!! coefficients \f$\alpha_P\f$ and \f$\alpha_I\f$ have dimensions
!! \f$t^{-1}\f$ and \f$t^{-2}\f$, respectively, and should be chosen
!! according to the dynamical characteristics of the problem. The
!! constant \f$\bar{\underline{M}}\f$ is defined in this module as \c
!! bar_mass_flow, and represents the desired mass flow.
!!
!! \todo The prescibed mass flow \f$\bar{\underline{M}}\f$ is now an
!! independent parameter; however, it is clear that the mass flow is
!! related to the mean velocity, and hence to the Reynolds number. So
!! the "free" parameters provided in this module are not really
!! independent. Using these relations, it should be possible to reduce
!! the number of independent parameters in this module.
!<----------------------------------------------------------------------
module mod_turb_channel_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_physical_constants, only: &
   mod_physical_constants_initialized, &
   t_phc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_turb_channel_test_constructor, &
   mod_turb_channel_test_destructor,  &
   mod_turb_channel_test_initialized, &
   test_name,        &
   test_description, &
   phc,              &
   ntrcs,            &
   bar_mass_flow,    &
   ref_prof,         &
   coeff_dir,        &
   coeff_visc,       &
   coeff_init,       &
   coeff_outref,     &
   coeff_f

 private

!-----------------------------------------------------------------------

 ! public members
 character(len=*), parameter   :: test_name = "turb_channel"
 character(len=100), protected :: test_description(3)
 type(t_phc), save, protected  :: phc
 integer, parameter            :: ntrcs = 0
 real(wp), parameter :: bar_mass_flow(3) = (/ 0.8_wp ,0.0_wp,0.0_wp /)
 character(len=*), parameter   :: ref_prof = "none"

 logical, protected :: mod_turb_channel_test_initialized = .false.

 ! private members
 real(wp), parameter :: &
   alpha = 0.7_wp,  & !< power low exponent for \f$\mu\f$
   pr    = 0.72_wp, & !< \f$Pr\f$
   gamma = 1.4_wp,  & !< \f$\gamma\f$
   ! coefficients of the forcing term
   alpha_p = 0.1_wp, &
   alpha_i = 0.5_wp
   ! the following variables are introduced for clarity
 real(wp), parameter :: &
   rho_r = 1.0_wp, &  !< \f$\rho_r\f$
   l_r   = 1.0_wp, &  !< \f$L_r\f$
   v_r   = 1.0_wp, &  !< \f$V_r\f$
   t_r   = 1.0_wp, &  !< \f$T_r\f$
   pi = 3.1415926535897932384626433832795029_wp, &
   box(3,2) = reshape(                                          &
     (/   0.0_wp  , -1.0_wp ,     0.0_wp       ,         &
        4.0_wp*pi ,  1.0_wp , 4.0_wp/3.0_wp*pi /),(/3,2/))
 real(wp) :: &
   re,     & !< \f$Re\f$
   ma,     & !< \f$Ma\f$
   re_tau, & !< \f$Re\f$
   delta0    !< amplitude of the initial random perturbation

 ! Input file ----------
 character(len=*), parameter :: &
   test_input_file_name = 'turbulent_channel.in'

 character(len=*), parameter :: &
   this_mod_name = 'mod_turb_channel_test'

 ! Input namelist, to be read from turbulent_channel.in
 namelist /input/ &
   re, ma, re_tau, delta0

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_turb_channel_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized          .eqv..false.) .or. &
       (mod_kinds_initialized             .eqv..false.) .or. &
       (mod_fu_manager_initialized        .eqv..false.) .or. &
       (mod_physical_constants_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_turb_channel_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'turb_channel, dimension d = ',d
   write(test_description(2),'(a)') &  
     'viscosity given by power law '
   write(test_description(3),'(a)') &  
     'input velocity: incompressible parabolic profile'

   !------------------------
   ! Redefinition of phc constants

   ! This test case does not consider gravity nor Earth rotation
   phc%omega   = 0.0_wp
   phc%gravity = 0.0_wp
   ! Air as an ideal gas
   phc%rgas    = 1.0_wp/(gamma*ma**2)
   phc%cp      = 1.0_wp/((gamma-1.0_wp)*ma**2)
   phc%cv      = phc%cp - phc%rgas
   phc%gamma   = gamma
   phc%kappa   = phc%rgas/phc%cp
   ! Sea level reference pressure (from the state equation)
   phc%p_s     = phc%rgas
   phc%p_sf    = 1.0_wp/phc%p_s
   ! Reynolds number
   phc%re      = re
   phc%re_tau  = re_tau  
 !------------------------

   mod_turb_channel_test_initialized = .true.
 end subroutine mod_turb_channel_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_turb_channel_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_turb_channel_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_turb_channel_test_initialized = .false.
 end subroutine mod_turb_channel_test_destructor

!-----------------------------------------------------------------------

 !> Kinematic viscosity:
 !! \f$\nu=\frac{\mu_0}{\rho}T^\alpha=\frac{T^\alpha}{\rho\,Re}\f$
 pure function coeff_visc(x,rho,p,t) result(nu)
  real(wp), intent(in) :: x(:,:), rho(:), p(:), t(:)
  real(wp) :: nu(2,size(x,2))

  nu(1,:) = (t**alpha)/(rho*re)
  nu(2,:) = (1.0_wp/pr)*nu(1,:)

 end function coeff_visc

!-----------------------------------------------------------------------

 !> Perturbed logarithmic profile
 !!
 !! The perturbation is added to the velocity profile, while keeping
 !! the reference values for density and temperature.
 !!
 !! To reduce the divergence of the perturbed flow, we have \f$\delta
 !! U_{x_1}=\delta U_{x_1}(x_2)\f$, \f$\delta U_{x_2}=\delta
 !! U_{x_2}(x_3)\f$ and so on.
 pure function coeff_init(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1)+ntrcs,size(x,2))

  integer :: id, idp
  real(wp) :: xm, dx, xi(size(x,2))

   ! density: uniform profile
   u0(1,:) = rho_r
   ! velocity: log profile with perturbations
   do id=1,size(x,1)
     idp = mod(id,size(x,1))+1
     xm = 0.5_wp*sum(box(idp,:))
     dx = box(idp,2)-box(idp,1)
     xi=(x(idp,:)-xm)*(2.0_wp/dx)
     u0(2+id,:) = delta0 * logistic_chaos(xi)
   enddo
   xi = (x(2,:)-0.5_wp*sum(box(2,:)))*(2.0_wp/(box(2,2)-box(2,1)))
   u0(3,:) = u0(3,:) + uu_par(xi)
   ! energy
   u0(2,:) = u0(1,:)*phc%cv*t_r + 0.5_wp*sum(u0(3:,:)**2,1)/u0(1,:)

   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e
 
 end function coeff_init

!-----------------------------------------------------------------------

 !> Velocity logarithmic profile for \f$y\in[-1\,,\,1]\f$
 !!
 !! We compute
 !! \f{displaymath}{
 !!  u^+ = \left\{\begin{array}{ll}
 !!   y^+, & y^+<5 \\[2mm]
 !!   \displaystyle \frac{1}{k}\log(y^+) + A, & y^+\geq 5 
 !!  \end{array}\right.
 !! \f}
 !! where \f$y^+=(1-|y|)Re_\tau\f$ and \f$u^+=\frac{Re}{Re_\tau}u\f$.
 !! The two values \f$k\f$, \f$A\f$ are chosen to ensure a smooth
 !! transition at \f$y^+=5\f$.
 elemental function uu_log(y) result(uu)
  real(wp), intent(in) :: y
  real(wp) :: uu

  ! local variables
  real(wp), parameter :: &
    k  = 0.41_wp, &
    a  = 5.2_wp
  real(wp) :: yp, up

   ! 1. wall unit y
   yp = (1.0_wp-abs(y))*re_tau
   ! 2. log and linear velocity profile in wall units
   if(yp.gt.5.0_wp) then
     up = log10(yp)/k + a
   else ! the profile is modified so that it is continuous for yp=5
     up = (log10(5.0_wp)/k+a)/5.0_wp * yp
   endif
   ! non-dimensional with respect to Re
   uu = re_tau/re * up

 end function uu_log

!-----------------------------------------------------------------------


!-----------------------------------------------------------------------

 !> Velocity parabolic profile for \f$y\in[-1\,,\,1]\f$

 elemental function uu_par(y) result(uu)
  real(wp), intent(in) :: y
  real(wp) :: uu

  ! parabolic laminar profile
  
     uu = 1.5_wp*(1.0_wp-y**2)
       
 end function uu_par

!-----------------------------------------------------------------------

 pure function coeff_outref(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1)+ntrcs,size(x,2))

   u0 = 0.0_wp

 end function coeff_outref

!-----------------------------------------------------------------------

 !> Most of the boundary condition is taken from the reference
 !! profile. However, we want to let the density fluctuate while
 !! keeping a constant temperature.
 pure function coeff_dir(x,u_r,u_i,breg) result(ub)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in), optional :: u_r(:,:), u_i(:,:)
  integer, intent(in), optional :: breg
  real(wp) :: ub(2+size(x,1)+ntrcs,size(x,2))

   ub(1,:)  = rho_r
   ub(2,:)  = rho_r*phc%cv*t_r
   ub(3:,:) = 0.0_wp

   if(present(u_i)) then ! fix T, leave generic rho
     ub(1,:) = u_i(1,:)
     ub(2,:) = ub(1,:)*phc%cv*t_r
   endif

   ! subtract the reference state
   ub(1,:) = ub(1,:) - u_r(1,:) ! rho
   ub(2,:) = ub(2,:) - u_r(2,:) ! e

 end function coeff_dir

!-----------------------------------------------------------------------

 !> Define the non-dimensional driving force vector.
 pure function coeff_f(x,mass_flow,imass_flow) result(f)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in) :: mass_flow(:), imass_flow(:)
  real(wp) :: f(size(x,1),size(x,2))

   ! Driving force acts only in the x direction
   f(1,:)  = (re_tau/re)**2                           &
            - alpha_p*(mass_flow(1)-bar_mass_flow(1)) &
            - alpha_i*imass_flow(1)
   f(2:,:) = 0.0_wp

 end function coeff_f
 
!-----------------------------------------------------------------------

 !> Generate a chaotic field \f$\eta(\xi)\in[-1\,,\,1]\f$, with
 !! \f$\xi\in[-1\,,\,1]\f$.
 elemental function logistic_chaos(xi) result(eta)
  real(wp), intent(in) :: xi
  real(wp) :: eta

  integer,  parameter :: n = 20       !< number of iterations
  real(wp), parameter :: r = 3.999_wp !< must be < 4
  integer :: i

   ! eta is rescaled on [0,0.5] and the following iterations will
   ! generate values in [0,1].
   eta = 0.25_wp*(xi+1.0_wp)
   do i=1,n
     eta = r*eta*(1.0_wp-eta)
   enddo
   ! rescale eta on [-1,1]
   eta = 2.0_wp*(eta-0.5_wp)
 end function logistic_chaos

!-----------------------------------------------------------------------

end module mod_turb_channel_test

